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ABSTRACT 

Stability analysis of thermocapillary convection in rectangular cavities is performed using direct numer- 
ical simulations. Influence of the Reynolds number(Re), the fluid Prandtl number(Pr) and cavity aspect 
ratio(Ar) on the motion is investigated. Neutral stability curves for transition to time- dependent con- 
vection are delineated in the Re - Ar plane for fluids with Pr=1.0, 4.4, 6.78 and 10. Several interesting 
features of these diagrams are discussed. One important conclusion is that Ar OT increases as Pr decreases. 
Thus, large values of both Ar and Re are necessary to induce thermocapillary oscillations for small Pr 
fluids such as liquid metals and semiconductor melts. Energy analysis is also performed for the oscillatory 
flow in the neighborhood of critical points in order to gain insight into the mechanisms leading to instability. 


INTRODUCTION 

Understanding fluid motion is crucial in some material processing technologies. In crystal growth from 
the melt, single crystals with uniform material properties are desired, but homogeneity in crystals can be 
destroyed if melt motion is unsteady [1]. In the terrestrial environment, buoyancy and thermocapillarity 
are two major causes for convection. However, in low gravity environment, thermocapillary convection 
becomes dominant [2]. 

Numerous experiments (for example, [3], [4] and [5]) have demonstrated the existence of instability 
of thermocapillary convection, i.e. when the Marangoni number(Ma) exceeds a critical value, the motion 
undergoes a transition from steady to oscillatory. 

Thermocapillary flows have received considerable interest. A rich body of numerical investigations 
are available in the literature (see [6], [7] and [8] ). Results of direct numerical simulation of oscillatory 
thermocapillary convection was reported in [9] by Peltier & Biringen. They provided a stability diagram in 
the (Ar, Ma) space for a Pr=6.78 fluid, and found a minimum critical Ar near 2.3 and a minimum critical 
Ma near 20,000 within the parameter range of Ar<3.8. 

Discussions of instability mechanisms can be found in [10] and [11] for dynamic thermocapillary infinite 
liquid layers, and in [12] for thermocapillary liquid bridges. Description of the oscillatory instability is also 
provided in [9], relating the temporal evolution of large-scale structures in the flow and their interaction 
with the temperature sensitive free surface. 

Here, we present a detailed stability diagram for fluids with Pr=10.0, 6.78, 4.4 and 1.0. Interesting 
features in the diagram are discussed. Comparison of flow patterns is provided to investigate the influence 
of Re and Ar on the motion of a Pr=10 fluid. In addition, energy analysis results are also given for 
convection of a Pr— 4.4 fluid with Re near both higher and lower critical points of the unstable region at 
Ar=3.0. 
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MODEL DESCRIPTION AND NUMERICAL PROCEDURE 


The physical model considered is thermocapillary convection of incompressible and Newtonian fluid in a 
rectangular cavity with height H and width At x H (At is the aspect ratio). Two vertical isothermal side 
walls are kept at T ^ on the left and T c on the right, respectively. Bottom boundary is rigid and adiabatic. 
Top boundary is a flat free surface open to a passive gas. Here, surface tension on the free surface is 
assumed to be a linear function of temperature as a — <tq - 7 (T - To). 

A dimensionless mathematical model in the stream function - vorticity formulation is used for numer- 
ical simulation, in which length, temperature, velocity and time are made dimensionless by use of scales 
H, AT = (T/j - T c ), 7 A T hi and H 2 /*/, respectively. Dimensionless parameters eu*e defined as: Pr = - 
and Re = 7 ~~ > where /x, v and a are dynamic viscosity, kinematic viscosity and thermal diffusivity, 
respectively. 

The coupled equation system is solved by a finite volume based scheme, in which the Poisson equation 
for the stream function is solved by the SOR method. Both the vorticity transport and energy equations 
are solved by the alternating direction implicit(ADI) method. All time derivatives and spatial derivatives 
including boundary conditions are approximated in second order accuracy. Velocities are obtained as spa- 
tial derivatives of stream function. Uniform mesh is used in the solution procedure with mesh resolution 
of 50 to 90 points per dimensionless unit length, depending on the Reynolds number considered. 

RESULTS 

For thermocapillary convection of a Pr— 6.78 fluid in rectangular cavities, Peltier & Biringen [9] con- 
structed a stability diagram in the (Ar, Ma) plane for the region Ar <3.8 and Ma<1.0xl0 5 . Their Ma is 
equivalent to our RePrAr. Some interesting characteristics were found, including the existence of double 
valued stability limits, i.e. as Ma goes up, the flow first changes from stable to oscillating at Ma^x , and then 
becomes stable again when Ma cr2 is reached. Ma cr2 grows monotonically with Ar, however, Ma^j does not. 

In the present work, we extend this investigation to fluids with Pr-1.0, 4.4 and 10.0, and construct 
the stability diagrams in the (Ar, Re) plane. A wider range of parameter space (0.0 < Ar < 7.0 and 0.0< 
Re < 1.3xl0 4 ) is considered as shown by Fig, 1, in which more interesting features are found. If we look 
at the particular fluid with Pr=4.4, the first critical aspect ratio is around 2.6. Unstable region exists for 
any Ar>2.6, and more interestingly, there are more than one unstable regions with Ar>6.0, i.e. if Re goes 
up from zero, one can find that the flow is steady at low Re, starts to oscillate at first critical point, goes 
back to steady state at second critical point, and becomes oscillatory again as Re reaches its third critical 
point. In addition, stability curves of fluids with different Pr do not cross each other. Neutral curves of 
smaller Pr fluids always locate inside curves of larger Pr, i.e. when Pr goes smaller, the critical aspect 
ratio always becomes larger, so does the lowest critical Reynolds number. From the trend given by these 
curves, we can draw a very important conclusion that, for fluids with very low Pr, large values of critical 
Ar and Re are expected for the transition to oscillatory thermocapillary convection. 

The convective flow field is strongly influenced by Re, Ar and Pr. For a Pr=10 fluid, Fig. 2 gives three 
examples of streamlines at Ar=3.0 and Ar— 6.0. (a) shows the flow field at Re^x for Ar=3.0, in which one 
can see a bi-cellular structure with a stronger cell near the hot wall and a much weaker cell close to the 
cold wall. Increasing Re to Re^ at about 7400, one can find, in (b), that the previous strong hot wall cell 
moves to the center of the cavity, and the weak cell disappears. For a larger aspect ration (Ar=6.0), Fig. 
2 (c) exhibits the flow pattern at Re cr i, from which we find that three cells exist, with the strongest one 
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still near the hot wall. Further increase of Ar will result in more cellular structure in the flow field. 

For convection of a Pr— 10 fluid at a large aspect ratio (Ar— 20), Fig. 3 displays the mean velocity 
profiles and snapshots of temperature fluctuation fields for three different values of Re— 1012, 1025 and 
1500, respectively. Comparison of these three mean velocity profiles shows almost identical patterns even 
when Re changes from 1012 to 1500, with most strong activities locating near two side walls. However, 
large difference can be found among the temperature fluctuation fields. At Re=1012 (b-1), which is very 
close to Re^i, a thermal wave generates near the center of the cavity and starts to die at the right edge of 
the strong flow cell near the hot wall. Most area in the right part of the cavity remains pretty calm. The 
wave actually propagates toward the hot wall if we look at snapshots at different time instants, which is in 
agreement with [10]. As Re goes up a little bit to 1025, (b-2) shows that the starting point of the thermal 
wave moves toward the cold wall. Further increase Re, this starting point keeps moving to the right until 
it reaches the cold wall. Fig. (b-3) shows the case when Re=1500, in which one can see a wave generating 
at the cold wall, propagating actively, and dying at the right edge of the strong flow cell near the hot wall. 


ENERGY ANALYSIS 

Energy analyses are performed for flows with Reynolds numbers in the neighborhoods of different criti- 
cal points. The physical parameters of the oscillatory flows are decomposed into their mean and fluctuating 
components, and investigations are conducted on the behavior of the fluctuation kinetic energy ( Ar) and 
the fluctuation thermal energy (6 — t 2 j 2). 

For the case Ar=3.0 and Pr=4.4, results of the energy analysis are provided here for flows with Reynolds 
numbers near both lower(Re=1950) and higher(Re=5020) critical points of the unstable region. Tempo- 
ral variations over one flow oscillation period are shown in Fig. 4(a-l) and (a-2) for the rate of change 
of the total fluctuation kinetic energy(dK jdr - ^(/ fcdfl), as well as its components 1^ (production), 
Ifc 2 (diffusion) and Ijt 3 (dissipation). It is seen that dK /dr oscillates with its time average being equal to 
zero, which means no kinetic energy is added to the flow over each period of oscillation. This is consistent 
with the fact that the flow field oscillates with a stable amplitude. If we look at I , Ifc 2 and I^ 3 , we 
find I and Ifc 3 providing two major contributions, with I always being positive (destabilizing) and U 3 
always being negative (stabilizing). Time averages of 1^, and I^ 3 are much larger than that of I^ 2 , however, 
the phase difference between I*., and I^ is always near 7r, which means that I^ 3 always cancels the effect 
of Ifc, . This gives the smaller term Ijt 2 a chance to influence the temporal behavior of dK /dr. In Fig. 
4(a-2), which is for the higher critical point, one can clearly see that dK/ dr oscillates at a very close am- 
plitude and a very small phase difference with I^ 2 , while the phase difference between Ifc, and I*. 3 is almost n. 

Variation of the rate of change of the total fluctuation thermal energy ( dQ/dr — j-(/ 6dVt) and its 
components (I f , and I t2 ) are given in Fig. 4(b-l) and (b-2) for the lower and higher critical points, 
respectively. As expected, dQ/dr oscillates with its time average being equal to zero, since the temperature 
field oscillates in a limit cycle with stable amplitude. In addition, although the time averages of l tl and 
It 2 have the same absolute value, the oscillation amplitude of I tl is much larger than that of l t2 . Thus 
It, dominates the oscillatory behavior of dQ/dr. In both cases, the phase difference between If, and 
dQ/dr is very small. If we further compare the magnitude of the destabilizing fluctuation thermal energy 
component (It, ) and the kinetic energy components (I*., and I fc 2 ), we find that the magnitude of thermal 
energy production I t| is generally two to three orders larger than that of Ijt, or I* 2 . Thus, I t , appears to 
be the major driving source of flow instability. 
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CONCLUSIONS 


Direct numerical simulation is employed for stability analyses of thermocapillary driven convection in 
rectangular cavities. Stability boundaries are delineated in the Re - Ar plane, in which several interesting 
features are found. Influence of Re and Ar on the flow patterns and the temperature fields is briefly 
discussed. In addition, energy analyses are performed to gain insight into mechanisms involved in the 
onset of instability. Results are presented for a Pr=4.4 fluid with Re near both higher and lower critical 
points at Ar=3.0. Investigations on 2D and 3D convection of lower Pr as well as larger Ar and Re are in 
progress. 
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Figure 2: Examples of streamlines for convection of a Pr=10 fluid. Case (a): Ar=3.0, Re=1530 (« the first 
critical number); Case (b): Ar=3.0, Re=7400 (~ the second critical number); Case (c): Ar=6.0, Re = 700 
(ss the first critical number). 
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Figure 3: Comparison of mean velocity profiles and temperature fluctuation fields of a Pr— 10 fluid at large 
aspect ratio Ar=20. Case 1: Re=1012, which is approximately the first critical point; Case 2: Re- 1025, 
which is slightly higher than the first critical point; Case 3: Re—1500, which is much higher than the first 
critical point. 
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Figure 4: Temporal variation of the rate of change of the total kinetic energy the thermal energy ( ~j~ 
and their components I fc] , I* 2 , U ;> , I £ , and I £ , for a Pr = 4.4 fluid at Ar = 3.0. ( a- 1 ) and (b-1) give results 
for Re=1950 (near the lower critical point), (a-2) and (b-2) provide results for Re = 5020 (near the higher 
critical point). 
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